home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libfft / TRY / r_check1du.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  4.9 KB  |  237 lines
  1. #include <stdio.h>
  2. #include <sys/time.h>
  3.  
  4. #include "fft.h"
  5. #include "constant.h"
  6.  
  7. /*                        */
  8. /*    Precision dependant declarations    */
  9. /*                        */
  10. #ifdef    DOUBLE
  11.  
  12. #define TOLERANCE   DTOLERANCE
  13. typedef        double    this_type;
  14.  
  15. #define        THIS_SUM    dsum_
  16. #define        THIS_GEN    dgen_
  17. #define        THIS_FTU    dft1du_
  18.  
  19. #define        THIS_FFTI    dfft1di
  20. #define        THIS_FFTU    dfft1du
  21. #define        THIS_SCAL    dscal1d
  22.  
  23. #endif
  24. #ifdef    SINGLE
  25.  
  26. typedef        float        this_type;
  27.  
  28. #define TOLERANCE   STOLERANCE
  29.  
  30. #define        THIS_SUM    ssum_
  31. #define        THIS_GEN    sgen_
  32. #define        THIS_FTU    sft1du_
  33.  
  34. #define        THIS_FFTI    sfft1di
  35. #define        THIS_FFTU    sfft1du
  36. #define        THIS_SCAL    sscal1d
  37.  
  38. #endif
  39.  
  40. /*                        */
  41. this_type THIS_SUM();
  42.  
  43. void inimat_();
  44. void GetArguments();
  45. void get_values();
  46.  
  47. int size, inc, n_trials;
  48. this_type *pa, *pb, *pRef, *pWork;
  49.  
  50. main(argc,argv)
  51. int argc;
  52. char *argv[];
  53. {
  54.     int i, ia, j, k, n_total, is_wrong, nerr;
  55.     double x, y, dx, dy, emax, sum, err, maxerr;
  56.     this_type    t;
  57.  
  58. /* ******************************************************* */
  59.     GetArguments( argc, argv);
  60. /* ******************************************************* */
  61.  
  62.     srandom( (123*getpid()) | 0x01);
  63.  
  64.     for( k = 0; k < n_trials ; k++) { 
  65.     get_values();
  66.     inc = 1;
  67.  
  68.     if( !(k%1) && (k)) printf("\n");
  69.     printf( "<%3d,%3d>", size, inc);
  70.     fflush(stdout);
  71.  
  72.     n_total = ((size+2)/2)*2;
  73.     pa = (this_type *)malloc( (3 * (n_total+2) + size + FACTOR_SPACE) * sizeof(this_type));
  74.     if( !(pa) ) {
  75.         fprintf( stdout, "Could not allocate ... Exiting\n");
  76.         exit (-1);
  77.     }
  78.     pb = pa + n_total + 2;
  79.     pRef = pb + n_total + 2;
  80.     pWork = pRef + n_total + 2;
  81.  
  82. /* *******************************************************
  83.     Initialize arrays
  84. ******************************************************* */
  85.     i = n_total+2;
  86.     THIS_GEN(pRef, &n_total);
  87.     pRef[n_total] = HUGE;
  88.     pRef[n_total+1] = -HUGE;
  89.     bcopy( pRef, pa, i*sizeof(this_type));
  90.     bcopy( pRef, pb, i*sizeof(this_type));
  91.  
  92. /*---NO SO PACKED VERSION---*/
  93. /* *******************************************************
  94.     direct Fourier Transform
  95. ******************************************************* */
  96.         j = -1;
  97.     THIS_FTU( &j, &size, pb, &inc);
  98.     pWork = THIS_FFTI( size, pWork);
  99.     is_wrong = THIS_FFTU( -1, size, pa, inc, pWork);
  100.  
  101.     emax = TOLERANCE*TOLERANCE*size;
  102.     err = 0.0;
  103.     sum = 0.0;
  104.     maxerr = 0.0;
  105.     j = 0;
  106.     for(i = 0, nerr=0 ; i < n_total ; i ++) {
  107.         t = pa[i] * pa[i];
  108.         x = pa[i] - pb[i]; x = x * x;
  109.         sum += t;
  110.         err += x;
  111.         if( x > maxerr) {
  112.         maxerr = x;
  113.         j = i;
  114.         }
  115.     }
  116.     err = sqrt(err / (double)size);
  117.     sum = sqrt(sum / (double)size);
  118.     maxerr = sqrt(maxerr);
  119.     printf(": avg:(%8.3g /%8.3g)= %8.3g <> maxerr: (%8.3g /%8.3g)= %8.3g \n\t", 
  120.     err, sum, err/sum, maxerr, sum, maxerr/sum);
  121.  
  122.     x =  THIS_SUM( &size, pRef, &inc);
  123.     dx = x - pa[0];
  124.     if( (t= dx *dx ) > TOLERANCE){
  125.         fprintf( stdout, 
  126.         "\nError direct FFT(%d) at index #0 : (%f)<->(%f)error(%f)",
  127.         size, pa[0], x, dx);
  128.  
  129.     }
  130.     dx = x - pb[0];
  131.     if( (t= dx *dx) > TOLERANCE){
  132.         fprintf( stdout, 
  133.         "\nError direct DFT(%d) at index #0 : (%f)<->(%f)error(%f)",
  134.         size, pa[0], x, dx);
  135.     }
  136. /* *******************************************************
  137.     inverse Fourier Transform
  138. ******************************************************* */
  139.     fflush( stdout);
  140.     is_wrong = THIS_FFTU( 1, size, pa, inc, pWork);
  141.  
  142.     i = 2 * inc;
  143.     j = (size -1) / 2;
  144.     x =  pb[0] + 2. * THIS_SUM( &j, pb+2*inc, &i);
  145.     if ( !(size & 1))
  146.         x += *(pb + (size) * inc);
  147.     dx = x - pa[0];
  148.     if( (t= dx *dx ) > TOLERANCE){
  149.         fprintf( stderr, "#");
  150.     } else {
  151.         fprintf( stderr, ".");
  152.     }
  153.  
  154.     t = 1. / (double)size;
  155.     THIS_SCAL( size, t, pa, inc);
  156.  
  157.     emax = TOLERANCE;
  158.     emax = emax * emax;
  159.     sum = 0.;
  160.     err= 0.;
  161.     nerr= 0;
  162.     maxerr= 0.;
  163.  
  164.     for(i = 0 ; i < (size*inc) ; i += inc) {
  165.         sum += (pRef[i] * pRef[i]) + (pRef[i] * pRef[i]);
  166.         x = pa[i] - pRef[i];
  167.         if( (t= x*x) > (emax))
  168.         nerr++;
  169.         err += t;
  170.         if( t > maxerr) maxerr = t;
  171.     }
  172.     err = sqrt(err / (double)size);
  173.     sum = sqrt(sum / (double)size);
  174.     maxerr = sqrt(maxerr);
  175.     printf(": avg:(%8.3g /%8.3g)= %8.3g <> maxerr: (%8.3g /%8.3g)= %8.3g", 
  176.     err, sum, err/sum, maxerr, sum, maxerr/sum);
  177.     if(nerr){
  178.         printf("\n%d errors detected \n");
  179.         exit(-2);
  180.     }
  181.     free(pa);
  182.     }
  183.     printf("\n");
  184.     return(0);
  185. }
  186.  
  187. int is_random;
  188.  
  189. void GetArguments( argc, argv)
  190. int argc;
  191. char *argv[];
  192. {
  193.     int i, j, k;
  194.     int nerror = 0;
  195.  
  196. #define ON    1
  197.  
  198.     n_trials = DEF_TRIALS;
  199.     is_random = 1;
  200.  
  201. /* ******************************************************* */
  202.     for ( i = 1 ; (i < argc) && (nerror != ON) ; i ++ ) {
  203.     if( argv[i][0] == '-') {
  204.         switch ( argv[i][1]) {
  205.         case 'i' :
  206.         case 'I' :  
  207.                 is_random = 0;
  208.                 break;
  209.         default  :  nerror = ON;
  210.         }
  211.     }
  212.     else {
  213.         if( i+1 > argc)
  214.         nerror = ON;
  215.         else { 
  216.         n_trials = atoi( argv[i]);
  217.         }
  218.     }
  219.     }
  220.     if( nerror == ON) {
  221.     fprintf( stderr, 
  222.         "Usage : %s [-p(arallel)] [n_trials]\n", argv[0]);
  223.     exit(-1);
  224.     }
  225. }
  226.  
  227. void get_values()
  228. {
  229.     inc = 1;
  230.   if( is_random){ 
  231.     size = random() % MAX_SIZE + 1;
  232.   } else{
  233.     printf( "Enter SIZE : ");
  234.     scanf("%d", &size);
  235.   }
  236. }
  237.